Spheromak formation and sustainment by tangential boundary flows 
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The nonlinear, resistive, 3D magnetohydrodynamic equations are solved numerically to demon- 
strate the possibility of forming and sustaining a spheromak by forcing tangential flows at the plasma 
boundary. The method can by explained in terms of helicity injection and differs from other helicity 
injection methods employed in the past. Several features which were also observed in previous dc 
helicity injection experiments are identified and discussed. 



Spheromak plasmas have been formed using several 
different techniques [TJ [2] . The existence of multiple for- 
mation methods clearly shows that the spheromak is a 
preferred (lowest energy) state toward which a magne- 
tohydrodynamic (MHD) system naturally evolves when 
the appropriate boundary conditions are imposed. The 
physical process responsible for the formation is mag- 
netic relaxation: on the time scale of MHD instabilities 
the plasma relaxes to the minimum energy state com- 
patible with its magnetic helicity content (which remains 
approximately constant) [3]. Once formed, the sphero- 
mak will decay in a resistive time scale because resis- 
tivity does dissipate magnetic helicity. For this reason, 
the sustainment of the configuration during times longer 
than the resistive decay time requires some helicity injec- 
tion method. Since relaxation operates on a shorter time 
scale, the spheromak configuration is maintained regard- 
less of the details of the specific helicity injection mech- 
anism. Some particular examples are the coaxial helicity 
injection method (CHI) |4-6 , the merging of helicity- 
carrying filaments (MHF) [7 and the helicity injected 
torus with steady inductive helicity injection (HIT-SI) 

In this Letter we report the first evidence coming from 
nonlinear, resistive, 3D MHD numerical simulations that 
demonstrate the possibility of forming and sustaining 
a spheromak by forcing tangential flows at the plasma 
boundary. The method can by explained in terms of he- 
licity injection and differs from other helicity injection 
methods employed in the past (CHI, MHF and HIT-SI). 

An enhanced helicity injection mode was recently re- 
ported in spheromak experiments with large plasma ro- 
tation [9 . Althought not analyzed in terms of boundary 
flows, this observation could support the feasibility of the 
mechanism proposed and studied in this Letter. 

We model the plasma using the resistive MHD equa- 
tions with finite viscosity and zero /3. The evolution equa- 
tions for u and B are: 

<9 £ u + uVu = (J x B)/p + • n, (1) 
d t B + V x E = 0, (2) 

where E=-uxB + r/J andn = ( Vu + Vu T ) - (2/3) x 
(V-u) (see Ref. [10] for further details). These equations 
are normalized with a (chamber radius), (imposed 



flux) and ca (Alfven speed). In addition, B and J are 
scaled with yT^o- Time is expressed in units of the Alfven 
time ta = cl/ca- The normalized resistivity r] and the 
kinematic viscosity v are set to 5 x 10 -5 . With these 
values the resistive time is r r ~ 800 (defined as in Ref. 
[TT]). The resulting system is solved with the Versatile 
Advection Code [12] . 

The domain is a cylinder of elongation h/a = 1, with 
perfectly conducting wall conditions (B n = and 
J x n = 0) and vanishing velocity (u = 0) at the cylin- 
drical boundary (r = a) and the top end (z = h). At 
the bottom end we impose the poloidal flux: ip(r,z = 
0) = C^ G {r/a) 2 (l - r/a) 3 , where C = 28.935 is a con- 
stant and ipQ is the maximum flux imposed by the gun. 
This geometry has been used to model the Spheromak 
Experiment (SPHEX) [13]. If u = is imposed at z = 
the full set of boundary conditions leads to vanishing he- 
licity and energy fluxes across the boundary. Recently, 
these conditions have been applied to study the decay of 
configurations representative of electrostatically driven 
spheromaks with open field lines [13] . 

Imposing tangential flows at a boundary intercepted 
by magnetic flux may result in the injection of helic- 
ity, as can be inferred from the equation dH/dt = 
-2j vV J- BdV-2<f dV [(A • B)(u • n)-(A • u)(B • n)]dS 
(neglecting electrostatic fields). The last term on the 
right gives the helicity injection produced by motions of 
the footpoints of the penetrating (open) magnetic field 
[T5] , Boundary shearing of magnetic fields has been ap- 
plied in the past to study reconnect ion events in low-/? 
plasmas relevant to the solar corona [T6] . 

The computations presented in this Letter have u r = 
u z = and uq = uq max[0, 25(2/5 — r)r], at z = 0. This 
flow produces helicity injection across the bottom end 
of the cylinder because Aq = ip(r,z = 0)/(2irr) there. 
The initial condition is the vacuum magnetic field and 
zero velocity inside the chamber. Setting uq = —0.1 we 
obtained the results shown in Figs. [TJ [2] and [3] The 
results of two more runs, with equal to —0.05 and 
—0.2, are discussed in the context of Figs. [4] and [5] 

The evolutions of the maximum poloidal flux (Vw) 
and the toroidal flux across the entire poloidal plane (<£) 
are shown in Fig. [I] (a). For t < 45, ijj ma remains at 
the constant imposed value, while there is a build up 
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FIG. 1: (Color online) Evolution of ip m a and <I> (a). In the 
inset, the evolution for t < 100 is showed in more detail. In 
addition the (halved) magnetic energy and the helicity are 
shown. Colormaps of A and ip contours in the poloidal plane 
are shown at three different times (b)-(d). 



of <£. Magnetic energy and helicity (computed using a 
gauge invariant definition [17 j) are shown in the inset of 
Fig. [l](a). At t ~ 45, a relaxation event which produces 
significant energy dissipation at approximately constant 
helicity and flux conversion from toroidal to poloidal is 
clearly observed. The kinetic energy K (not shown) is 
small during the whole evolution (K/W has a peak of 0.1 
at t ~ 50 and remains below 0.01 during sustainment). 
Panels (b), (c) and (d) of Fig. [I] show contours of i/j and 
colormaps of A (where A = J • B/£? 2 , is computed using 
the n = component of the toroidal Fourier decompo- 
sition of B) at different times. At t = 30, A is strongly 
peaked near the geometric axis and no closed ijj contours 
are identified. After the relaxation event (t = 70) A is re- 
distributed and closed ip contours appear, indicating the 
formation of a spheromak configuration. The poloidal 
and toroidal fluxes continue to increase until t ~ 700 and 
afterwards a quasi stationary state is sustained for one 
resistive time, indicating that a balance between helicity 
injection and dissipation has been reached. The ip and 
A distributions shown in Fig. [I] (d) are representative of 
this quasi-steady state. 

It is well known that spheromak formation and sus- 
tainment necessarily involves non-axisymmetric activity. 
The magnetic field lines (followed from fixed positions 
at the bottom end) plotted in Fig. [5] clearly show that 




FIG. 2: (Color online) Twenty field lines followed from fixed 
positions at the bottom end, for eight different times. The 
opacity of the lines is gradually decreased when their length 
L is longer than 4, and they become completely transparent 
when L > 5. The structure observed at t = 110 is almost in 
its final saturated state (observed at t = 1072) and rotates 
clockwise (when viewed from top). 



our computations reproduce this feature. The initial 
field lines of the vacuum solution (t = 0) expand and 
twist (t = 8) forming a central current-carrying column 
(t = 16). As the current through the central column in- 
creases and the field lines increase their twisting, the con- 
figuration eventually becomes unstable (t = 37) and the 
typical helical structure of a kink instability quickly de- 
velops (t = 46). This instability, with dominant toroidal 
wavenumber n = 1, saturates at a relatively small and 
fixed amplitude (W n=1 /W n= o ~ 0.01 during sustain- 
ment, see Fig. [4| causing the central column to adopt 
an almost fixed helical shape which rotates as indicated 
in Fig. [2] (t = 110 and t = 1072). Our results repro- 
duce the experimental observation of a fixed rotating he- 
lical structure with a strongly asymmetric return column 
[T8] . In contrast with our case (driven purely by bound- 
ary flows), the plasma in those experiments was driven 
electrostatically. We note, however, that the presence of 
an inductive component in the electric field when mag- 
netic fluctuations were active, was also reported [18] . 

It is important to note that the motion of this struc- 
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FIG. 3: (Color online) Twenty instantaneous profiles (be- 
tween ti and £3) of Be and B z , at z = h/2 and = 0. The 
profile of the time averaged fields is also shown. The n — 
component of the toroidal velocity at £2 is shown in the inset. 



ture is produced by the coherent oscillation of the n = 1 
mode and not by a rigid rotation of the plasma. This is 
shown in Fig. [3] where the radial profiles of Bq and B Z1 
at z — h/2 and = 0, are plotted for several times be- 
tween £1 = 1290 and £3 = 1380. The profiles at £1 and £3 
are almost identical, indicating that the structure makes 
one turn every 90 Alfven times. Such rotation is much 
slower than that expected from the imposed uo = —0.1 
at r = 0.2. Furthermore, the plasma toroidal velocity 
reverses its sign as can be observed in the inset of Fig. 
[3] (this velocity colomap is plotted at £2 = 1335 but it 
is representative of the flow pattern during sustainment, 
700 < £ < 1500). This velocity reversal is also observed 
in the other two runs presented in this work and it is ev- 
idently a feature of the saturated state of the instability. 
However, the specific reason for this velocity reversal is 
still not yet understood. 

The evolution of %jj ma and the poloidal magnetic field 
near the wall (total B z \ w and B™ =0 \ w ) at z = h/2 and 
= are shown for three values of uq in Fig. [4] (a)-(c). 
The magnetic energy (relative to the initial energy) of 
the modes n = 0, n = 1 and n = 2 is shown in Fig. [4] 
(d). These results indicate that increasing the helicity 
injection rate (which is linear in uq) leads to configura- 
tions with higher flux amplification and more energy in 
the n — mode. Nevertheless, the energy associated to 
the n > modes (in particular the n = 1) saturates at 
a fixed amplitude, which is roughly independent of uq. 
This implies a lower level of fluctuations relative to the 
n = mode at higher helicity injection rates. A similar 
saturation mechanism has also been observed in experi- 
ments [151115]. 
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FIG. 4: (Color online) Evolution of Vw and the poloidal 
magnetic field near the wall (total B z \ w and B^°\ w ) for three 
values of uo (a)-(c). Field lines showing the central open flux 
column are presented in each case at the time indicated by 
the arrow. Magnetic energy of the modes n — 0, 1 and 2 (d). 



The oscillation of the fluctuation in the poloidal mag- 
netic field near the wall (B z \ w ) is clearly observed in Figs. 
|4] (a)-(c). As discussed previously, this corresponds to a 
coherent oscillation which produces the rotation of the 
magnetic structures shown in the insets of Figs. [4] (a)- 
(c). Note that the rotation frequency increases rapidly 
with the helicity injection rate. Another important ob- 
servation is that the rotation frequency is not the same, 
in general, that the frequency of the oscillating dynamo 
term that amplifies ipm a (which depends on the correla- 
tion of the fluctuations of u and B). This is specially 
clear in Fig. [I] (b). 

The MHD activity maintains the configurations ob- 
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FIG. 5: (Color online) X(ip) profiles for the case with uq — 
—0.1 at ten times between t\ and £3 (a). The kink instability 
threshold is indicated. Time evolution of A at the geometric 
axis (Xga) and at the magnetic axis (Am a), for the three runs 

(b). 

tained for the three different values of uq close to the 
kink stability boundary. This is shown in Fig. |5j The 
A(^) profiles of the n = mode are plotted in Fig. [5] 
(a), for the case with uq = —0.1, at ten times between t\ 
and £3 (ip is normalized with Vw)- The indicated kink 
instability threshold is a linear A(^) profile, which has 
a slope a = —0.4 [20]. Fig. [5] (b) shows the evolution 
of the A values at the geometric axis (Xga) and at the 
magnetic axis (Xma) f° r the three runs. The behavior 
observed here is similar to that described in previous ex- 
periments [19, 20 . The configuration evolves around the 
kink stability boundary as a result of the competition be- 
tween the external forcing, which tends to steepen the A 
profile, and the relaxation, which tends to flatten it. 

In summary, we demonstrated the possibility of form- 
ing and sustaining a spheromak by imposing tangential 
flows at the boundary. Several aspects of the process were 
described. The relationship between our results and the 
existing experimental evidence was discussed. The fact 
that our simulations reproduce many features observed 
in electrostatic CHI experiments is interpreted as follows. 
Since the spheromak is formed and sustained by the re- 
laxation of an unstable configuration, the dynamics of 
the process should be independent of the details on how 
this configuration is driven unstable. In this sense, this 
Letter not only provides evidence on a new spheromak 
formation and sustainment mechanism, but it also pro- 



vides valuable information pertaining to the dynamics 
of the kink instability in electrostatically driven sphero- 
maks. To conclude, we note that althought our results 
were obtained for a spheromak, boundary plasma flows 
could also be used to inject helicity in other configura- 
tions (i.e. spherical tokamaks). 
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